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Abstract 

We consider a saturated porous medium in the solid-fluid segregation regime 
under the effect of an external pressure applied on the solid constituent. We 
prove that, depending on the dissipation mechanism, the dynamics is described 
either by a Cahn-Hilliard or by an AUen-Cahn-like equation. More precisely, 
when the dissipation is modeled via the Darcy law we find that, provided the 
solid deformation and the fluid density variations are small, the evolution equa- 
tion is very similar to the Cahn-Hilliard one. On the other hand, when only 
the Stokes dissipation term is considered, we find that the evolution is governed 
by an AUen-Cahn-like equation. We use this theory to describe the forma- 
tion of interfaces inside porous media. We consider a recently developed model 
proposed to study the solid-liquid segregation in consolidation and we give a 
complete description of the formation of an interface between the fluid-rich and 
the fluid-poor phase. 

Keywords: phase transformation (A), porous material (B), finite differences 
(C), energy methods (C) 



1. Introduction 

Many interesting swelling and shrinking phenomena are observed when a porous 
medium deforms as a consequence of the variation of the fluid mass content. 
The behavior of the system can be controlled via different external parame- 
ters, e.g. the fluid chemical potential [Gelb et al. 19991 IConvertino et al. 20031 
IConvertino et al. 2002] or a mechanical pressure exerted on the solid component 
[CiriUo et al. 20111 ICmllo et al. 2010] . 

By focusing the attention on this last class of phenomena we are led to con- 
sider the so-called consolidation problem. Despite the existence of classical and 
well known theories |von Terzaghi 1946[ IBiot 1941[ [Cryer 1963| IMandel 1953] . 
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the dilatant/contractant behavior of porous soUds needs, however, deeper in- 
vestigation. We refer, in particular, to undrained conditions, i.e., to the case 
in which the loading rate does not allow pore-water pressure dissipation. In 
these conditions, the coupling between shearing and dilatancy typically induces 
the formation of strain localization bands. In particular shear bands with tran- 
severse size depending on the texture and the constitutive properties of the 
material, see e.g. [Besuelle 2001j . are seen. For instance, in simple shear test, 
densely or loosely packed granular media exhibit dilatancy or compaction behav- 
ior. Moreover due to porosity change, the fluid can migrate through the pores 
and eventually remain segregated, possibly enhancing localised overpressuriza- 
tion and fluidization of the soil, see e.g. [Kolymbas 1994[ INichols et al. 1994] . 

Local dilatancy is typically modeled in the framework of (poro-) plasticity by 
assuming that irreversible deformations take place inside the porous medium, see 
e.g. the general treatises [Coussy 2004[|Coussy 2010] . Special important appli- 
cations to multi-phase porous media and geomaterials have been considered in 
the recent literature [Zhang et al. 1999[[Coriin et al. 2006llChambon et al. 20011 
lAndrade et al. 20TT] by approaching plasticity with strain gradient [Fleck et al. 19971 
IFleck et al. 2001] or multi-scale analysis [Christoffersen et al. 198l] |Ncmat Nasse r et al. 1993] . 

In the papers [Cirillo et al. 20111 ICirillo et al. 20091 ICirillo et al. 2010] we 
have attacked this problem from the point of view of bifurcation theory and 
we have shown that it is possible to describe interesting phenomena (still in 
the range of non-linear elasticity) taking place when the pressure exerted on 
the solid exceeds a suitable limiting value. More precisely, we have proven that 
enhanced models based on gradient elasticity allow for describing the transition 
between phases of the porous medium associated with different fluid content, 
say a fluid-rich and a fluid-poor phase. In other words the theory proposed in 
those papers is able to describe the occurrence of segregation and local overpres- 
surization of the pore-water, see the above quoted references for more details. 

Experiments demonstrate that the presence of a pressurized fluid-rich phase 
at the base of a layered granular material can give rise to an unstable (with re- 
spect to possible fluidization) conflguration of the system [Kolymbas 1998[ . This 
result is achieved in laboratory, see [Nichols et al. 1994] , considering a fluidized 
column test setup where a fluid is forced to flow through a saturated sample 
from the bottom. By tuning the velocity of the fluid, the drag force acting 
on the solid grains, and hence the possible unbalance of gravity, is controlled 
[Vardoulakis 2004a| IVardoulakis 2004b[ . The model in [Cirillo et al.^009] is 
able to explain the formation of fluid-rich regions inside a porous material via a 
sort of segregation process, due to different initial data and parametrized by the 
pressure applied to the solid; note that the control parameter in the experiment 
quoted above is, instead, the velocity of the injected fluid. Different dissipation 
mechanisms, driven by the microstructural properties of the material, give rise 
to different deformation paths and fluid mass content evolutions. 

In this paper we show that, depending on the dissipative mechanism which 
is taken into account, the evolution is described either by an AUen-Cahn or 
by a Cahn-Hilliard-like system of partial differential equations, which indeed 
play an important role in the theory of phase ordering dynamics [Bray 1994} 
[Langer 1992[ [Eyre 1993] . When a system is quenched from the homogeneous 
phase into a broken-symmetry one, think to a ferromagnet or to a gas abruptly 
cooled below their critical temperature, the two phases have to separate and 
order has to increase throughout the system via domain coarsening. This phe- 
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nomenon can be described via a field u on the physical space C M**, with d the 
physical dimension, having the proper physical interpretation, for instance local 
magnetization in fcrromagncts, density in liquid-vapor systems, and concentra- 
tion in alloys. A possible model for the evolution of this field is the Allen-Cahn 
equation 

(1) 

where e is a positive constant and W{u) a double well regular function. The 
Allen-Cahn equation, also called the time-dependent Ginzburg-Landau equa- 
tion, was introduced in [Allen et al. 1979] to describe the motion of anti-phase 
boundaries in crystalline solids. In this context u represents the concentration 
of one of the two components of the alloy and the parameter the interface 
width. 

The Allen-Cahn equation is not suited to describe phase separation when 
the order parameter is conserved, namely, when the integral of the field u on the 
whole space Vl is constant. On the other hand this is possible in the framework 
of the Cahn-Hilliard equation 

-^-^{e^/^u~W'{u)) (2) 

by assuming di^u = d^l^u = on d^l, that is by requiring the derivative of 
u and that of the Laplacian of u orthogonal to the boundary to vanish on 
the boundary dVl. itself. The Cahn-Hilliard equation was firstly introduced to 
describe spinodal decomposition in alloys |Cahn 19^T| . 

A straightforward way to derive the Allen-Cahn and the Cahn-Hilliard equa- 
tions is that of assuming the evolution of the field u to be governed by a gradient 
equation [Fife 20021 lAlikakos et al. 1991| du/dt = —SF/Su associated with the 
Landau energy functional 

F{u):= [ \le^\\Vu\\^ + W{u)\dx (3) 
JO '-2 J 

The functional above has a clear physical interpretation: the term W has two 
minima corresponding to the two phases, while the gradient-squared term asso- 
ciates an energy cost with an interface between the two phases. If no constraint 
to the total value of the order parameter is imposed, it is possible to compute 
the gradient of the Landau functional in the Hilbert space L'^{n) to get the 
Allen-Cahn equation. When the integral of the field u is assumed to be con- 
stant throughout the evolution, computing the gradient in the L'^{Q) results 
into a non-local evolution equation, while by using the Hilbert space H^^(fl) 
the Cahn-Hilliard equation is found. 

The above described approach is general and applies to any physical inter- 
pretation of the two equations. When the discussion is bounded to particular 
physical systems, phenomenological approaches can also be used requiring the a 
priori specification of the constitutive equations. The related physical literature 
is huge and, among others, we refer to [Gurtin 1996j . In the context of binary 
fluids, for example, u is the fluid concentration. Assume the continuity equation 
in the form + V • J = 0, where J is the fluid current prescribed in terms of 
the chemical potential by Pick's law J = — fcV/i, with k a positive constant. 
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By letting /i be the functional derivative of F, one gets 

— = _V • J = -V • (-fcV/i) = kAfi = yfcA— = kA(-e^Au + W'(u)) 
ot ou 

which is the Cahn-Hilliard equation. 

In our problem the AUen-Cahn and the Cahn-Hilliard-like equations ap- 
pear in a completely natural way by using the variational description of con- 
tinuum mechanics. By following [Coussy 2004[ ISciarra et al. 2008j we describe 
the porous material via two fields, the strain e and the fluid density (measured 
with respect to the solid reference volume) m. We then assume that the possi- 
ble motions of the system are those such that the variation of the action equals 
the opposite of the time integral of the virtual work of the dissipative forces 
corresponding to the variation of the fields. Solving the variational problem we 
get the equations of motion. 

The explicit form of the equations of motion depends not only on the form 
of the action but also on the way in which the dissipation phenomenon is mod- 
eled [Nield 20071 |Lopatnikov et al. 2004| . Assuming the potential energy to be 
quadratic in the first derivatives of the strain and of the fluid density, the evolu- 
tion is described by an Allen-Cahn-like equation for the field m provided that 
dissipative forces are assumed to be proportional to the second derivative of the 
velocity of the fluid. We refer to this assumption as to the Stokes dissipation. 
On the other hand, we flnd that the evolution is described by a Cahn-Hilliard- 
like equation for the field m provided dissipation is modeled via the Darcy law, 
namely, the dissipative forces are proportional to the velocity of the fluid. 

It is possible to give a nice physical interpretation to the fact that the evolu- 
tion is described by the AUen-Cahn or the Cahn-Hilliard-like equations depend- 
ing on the way in which dissipation is put into the game. With homogenization 
arguments it has been proven that when the characteristic length at the local 
{£) and at the macroscopic (L) scales separate, say £/L <^ 1, then the resis- 
tance experienced by the fluid when flowing through the porous material can 
be described by the Darcy law [Ene et al. 1975j . On the other hand when the 
separation of scales is poor, the Darcy limit unfairly approximates the flow, 
since the macroscopic characteristic length is of the same order of magnitude 
as the pore characteristic size [Auriault et al. 2005j . In other words in this case 
the solid grains have to be considered as microporous themselves. 

Moreover, we note that the steady Stokes equations, which describe the be- 
havior of the fluid at the microscopic level, can be upscaled to a set of equations 
describing a non-Stokesian flow only if the size of the solid obstacles to the fluid 
motion have a critical size with respect to the inter-obstacle distance, see for 
more details [Allaire 19971 Theorem 3.1]. This condition provide the so-called 
Brinkman limit. Conversely when the ostables are too small the homogenized 
equations coincide with the Stokes equations. Thus the Allen-Cahn-like equa- 
tions, which correspond to a purely Stokesian flow, are controlled by the vis- 
cosity of the fluid; on the other hand the Cahn-Hilliard-like equations, which 
correspond to a Darcy flow, are controlled by the permeability of the solid skele- 
ton. In particular the Cahn-Hilliard-like model reads, within the framework of 
soil consolidation theory, as a generalization to strain gradient materials of the 
one-dimensional Terzaghi problem, see |Madeo et al. 200"8] . 

The paper is organized as follows. We first introduce the model in Section [2] 
In Section [3] we assume small variations of the fields with respect to some ref- 
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erence values, and we find the Allen-Cahn and Cahn-Hilliard-like equations. 
In Section [4] we consider a special model allowing to describe solid-fluid segre- 
gation in consolidation jCirillo et al. 2011l[Cirillo et al. 2009llCirillo et al. 2010] 
and we give a numerical study of the equations deduced in the previous sections. 

2. The model 

In this section we introduce the one dimensional poromechanical model whose 
geometrically linearized version will be studied in the following sections. Kine- 
matics will be briefly resumed starting from the general statement of the model 
[Coussy 2004| together with some particular issue introduced in [Sciarra et al. 2008] . 
The equations governing the behavior of the porous system will then be deduced 
prescribing the conservative part of the constitutive law through a suitable po- 
tential energy density <i> and the dissipative contributions through purely Stokes 
or Darcy terms. Special emphasis will be given to the boundary conditions and 
the extended definition of the essential and the natural ones will be discussed. 

2.1. Poromechanics setup 

Let Bs := [^1,^2] C M, with £i,£2 E M, and := M be the reference configu- 
rations for the solid and fluid components [Coussy 2004| . The solid placement 
Xs : -Bs X M — ^ M is a function such that the map xs(-,t), associating to 
each Xs G -Bs the position occupied at time t by the particle labeled by Xg 
in the reference configuration Bs, is a C^-diffeomorphism. The fluid place- 
ment map : i3f X M — >■ R is defined analogously. The current configuration 
Bt ■— Xs{Bs, t) at time t is the set of positions of the superposed solid and fluid 
particles. 

Consider the function (p : Bg x M. B^ such that (f>{Xs,t) is the fluid 
particle that at time t occupies the same position of the solid particle Xs', 
assume, also, that 4>{-,t) is a C^-diffeomorphism mapping univocally a solid 
particle into a fluid one. The three fields xs, Xf j and </> are not at all independent; 
indeed, by definition, we immediately have that Xi{4>{^Si't)-,'t) — Xs{^Si't) for 
any Xs € Bs and t e R. 

In the sequel we shall often use the inverse functions of the space sections 
of the field Xf, XS: and cf). We shall misuse the notation and let (j)~^{-,t) be 
the inverse of the map Xs cl){Xs, t) at a given time t. Similarly we shall also 
consider Xs^{.-,t) and Xf ^'^O- 

The Lagrangian velocities are two maps associating with each time and each 
point in the solid and fluid reference space the velocities of the corresponding 
solid and fluid particles at the specified time. More precisely, the Lagrangian 
velocities are the two maps Ua : B^ x K — >■ K defined by setting 

u^{X^,t):=^{X„,t) (4) 

for any Xa G B^, where a — s,f. We also consider the Eulerian velocities 
Va : i?t X M -> E associating with each point x E Bt and for each time t eR the 
velocities of the solid and fluid particle occupying the place x at time t; more 
precisely we set Va{x,t) := WQ(xQ^(a;, t), t). 

In studying the dynamics of the porous system one can arbitrarily choose 
two among the three fields xsj Xfj and (/)■ Since the reference configuration Bs 
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of the solid component is know a priori, a good choice appears to be that of 
expressing all the dynamical observables in terms of the fields xs and cf) which 
are defined on Bs- Consider, for instance, the Lagrangian velocity of the 
fluid component which is defined on x M; we prove that for any Xg G Bg and 

^f(0(Xs,i),<) =xs(^s,i)- i^^Xs(^s,t) (5) 

where the dot and the prime denote, here and in the sequel, the partial derivative 
with respect to time and to the space variable Xg respectively. The above 
equation gives the expression of the fluid Lagrangian velocity in terms of xs 
and (f). 

In order to prove (5|, note that X{{^{i — Xs{4'~^{X{j t)i t) for any X^ and 
t, where we have used~the definition of fluid placement map. By the definition 
Q of the Lagrangian velocity of the fluid, we then get 

Wf(Xf, t) = Xs(0"' (% 0, or H^f, t) + Xs(0"'(^f, t),t) 

Since <f>^^{-,t) is the inverse function of ((>{■, t), we have that 4>{(t)^^{X-^,t),t) = 
X^ for any X^ G B^. By deriving this equality with respect to time we get 

(i>'{r\x^,t),t)(^-^ix^,t) + <j){r\Xi,t),t) = 

The last two equations yield ([s]). 
2.2. Variational principle 

In order to write the equations of motion of the system we use a variational ap- 
proach much similar to that developed in |Sciarra et al. 2008] ; the differences be- 
tween the two computations will be pointed out. For the sake of self-sufficiency 
we shall report the computation in detail for the one-dimensional case. 

It is natural to assume that, if the system is acted upon only by conservative 
forces, its dynamics is described by a Lagrangian density relative to the solid 
reference configuration space volume, depending on the space variable Xg and 
on time through (in principle) xs, 0, Xsi 4'"^ Xsj '^'j XSj and <p. The Lagrangian 
density is equal to the kinetic energy density minus the overall potential energy 
density accounting for both the internal and the external conservative forces. 
The equations of motion for the two fields xs and (p can be derived assuming 
that the possible motions of the system in an interval of time (ii,i2) C M are 
those such that the fields xs and (/) are extremals for the action functional 



A(xs,0):=/ dt dXs^(xs(^s, *),..., 0(^s,O) (6) 

in correspondence of the independent variations of the two fields xs and (j) on 
Bs X (ti,t2). In other words any possible motion of the system in the consid- 
ered interval is a solution of the Euler-Lagrange equations associated to the 
variational principle 6 A = 0. 

A different variational principle is needed if the fluid component of the system 
is acted upon by dissipative forces; the virtual work made by these forces must be 
taken into account. Consider the independent variations 5xs and Scf) of the two 
fields Xs and (f) and denote by 6W the corresponding elementary virtual work 
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made by the dissipative forces acting on the fluid component. The possible 
motions of the system, see for instance [Blanchard et al. 19921 Chapter 5] , in an 
interval of time (ti, ^2) C M are those such that the fields xs and cj) satisfies the 
variational principle 

SA^- SWdt (7) 

namely, the variation of the the action integral in correspondence of a possi- 
ble motion is equal to the integral over time of minus the virtual work of the 
dissipative forces corresponding to the considered variation of the fields. 

Up to now the discussion has been conducted along very general lines. In the 
following sections we shall specify first the form of the virtual work of the drag 
forces and then that of the action. We shall obtain, finally, an explicit repre- 



sentation of the equations of motion (see (26 1 and (27)) with suitable boundary 
conditions. 

2. 3. Virtual work 0/ the dissipation forces 

The way in which dissipation has to be introduced in saturated porous media 
models is still under debate, see for instance [Nield 2007) . In particular accord- 
ing to the effectiveness of the hypothesis of separation of scales, between the 
local and macroscopic level, Darcy's or Stokes' effects are accounted for. 

We first consider the so-called Darcy effect, i.e., the dissipation due to forces 
proportional to the velocity of the fluid component measured with respect to the 
solid. By following phenomenological arguments [Bear 1972] or by developing 
suitable homogenization schemes |Ene et al. 1975j it is seen that the fluid flow 
is, in this case, controlled by the permeability of the porous material. A natural 
expression for the virtual work of the dissipation forces acting on the fluid 
component and taking into account the Darcy effect is 

5Wj^ := - I D[v^{x,t) - vs{x,t)][5xi{x^\x,t),t) ~ 5xs{Xs\^,t).t)]dx (8) 

J Bt 

where _D > is a constant proportional to the inverse of permeability and 5x{ 
is the variation of the field X{ induced by the independent variations 5xs and 
54>. The quantity ~D[v^ — ws] is the dissipative force density (relative to the 
current configuration space volume) which depends on the kinematic fields only 
through the velocity of the fiuid component relative to the solid. Apparently 
equation ([s]) for the virtual work of Darcy dissipative forces is consistent with 
the classical expression of Darcy dissipation associated with the viscous flow 
through a porous continuum, see [Coussy 2004| . 

Following the recipe described above we have to express the virtual work in 
terms of the field xs and (j). We first remark ISciarra et al. 2008} Appendix C] 
that 

<5xf (0(^s, t),t)- SxsiXs. t) = -^^^ ^^(Xs, t) (9) 

for any Xs e i?s and i g E. Moreover, noted that usiXs^t) ~ xsiXs^t), by 
using the definition of the Eulerian velocities, from ^ we get 

v^{x,t) - vs{x,t) = (10) 



7 



In order to get an expression of the virtual work comparable to that of the action, 
we have to rewrite the integral as an integral extended to the solid reference 
configuration. Given t, by performing the change of variables x = Xsi^Sji) in 
the integral in ([s]) and by using equations ^ and ( 10 ) we get 



"■■^ ^ -"L litl»^<-''»-'> fill "»^<-''^-'>'^^ 

For the virtual work of the drag forces it is often considered also the contri- 
bution of the Stokes effect |Nield 2007] . namely, the virtual work of dissipative 
forces controlled by the second derivative of the velocity of the fluid component 
relative to the solid. This contribution is typically the leading term in the case 
when separation of scales between the local and the macroscopic level is no 
more valid. In these cases dissipative stresses must be taken into account at 
macroscopic level too [Allaire 1997j and, thus, the second derivative of the fluid 
velocity appears in the governing equations. A naive choice for the virtual work 
of the dissipation forces would be 

Siv^ix, t) - vs{x, t)]" [Sxi{xi\x, t),t)- Sxs{Xs\^, (12) 

Bt 

with S a positive constant; note that the second derivative is taken with re- 
spect to the space variable x in the current configuration Bt- Following the 
same thermodynamic argument as that developed for Darcy dissipation, to the 
expression above of the virtual work corresponds the effective power 

S[v^{x, t) - vs{x, t)]" [v^{x, t) - vs{x, t)] dx (13) 

Bt 

Since it is not negative definite, the proposed expression of the virtual work is not 
physically acceptable. This problem has been nicely solved in jSciarra et al. 2008j 
where the following expression of the virtual work 

<5%:=- / S[v^[x,t)-vs{x,t)n5xi{x'i\x,t),t)-5xs{Xs\x.t),t)]' (14) 



with > 0, has been proposed. Equation ( 14) for the virtual work of dissipative 
forces is indeed definitely consistent with the expression of the fluid dissipation 
when the fluid velocity satisfies the Navier-Stokes equations. As a matter of 
facts replacing the virtual displacement with velocities in the above expression 
yields the effective power of the Stokes drag forces which results to be negative 



definite. Note that, by performing an integration by part of ( 14 1, a term similar 



to (12 1 is found plus a contribution on the boundary of the actual domain. 



We rewrite, now, the expression of the Stokes virtual work as we have done 



for the Darcy contribution and we shall find an equation similar to (11). In 
other words we shall express all the functions appearing in ( [T4| in terms of the 
fields xs E^nd (p and, via the changing of variable x = xs(A's,i) for any fixed 
t, extend the integral to the domain B^. The computation is performed in the 
[Appendix A[ the result is 

'dXs 
(15) 



6Ws^-Sj^ — (#'x^' + 0>'xs - #"Xs) [S^' + ^ (Xs>' - 0"Xs)'5<^ 



8 



where we have omitted, for simpUcity, to write expUcitly that aU the functions 
are evaluated in {Xs,t). 



From (111 and ( |15[ ) we finaUy have the following expression for the virtual 
work of the drag forces obtained by taking into account both the Darcy and the 
Stokes effects 

5W := SWjy + SWg = - / {R5(j) + QS(f>') dXs (16) 



R Dj^^ix'sf + S^^.iWx'i + ^'^'X's " '^0"Xs)(Xs'^' " ^"x's) (17) 



where we have set 



and 

Q = ^^(00'xs + '/''0'Xs - '/'0"Xs) 
Since we need an expression of the virtual work in terms of the variations 



(18) 



5xs and 5(l> of the basic fields, we integrate ( 16 ) by parts and get 

5W = -{Q5^)%-[ {R~Q')S^dXs (19) 
Jbs 

As already noticed only in the Brinkman limit Darcy and Stokes dissipative 
effects can be accounted for together; on the other hand, in the general case, 
they are treated separately. This is indeed what we shall do in Section [3] 

2.4- Variation of the action 

In order to write explicitly the variation of the action we specify, now, the form of 
the Lagrangian density. In the sequel we shall not consider the inertial effects, so 
that, the Lagrangian density will be the opposite of the potential energy density 
associated to both the internal and external conservative forces. As it has 
been shown in [Sciarra et al. 2008] (see equation (18) therein) it is reasonable 
to assume that the potential energy density depends on the space and time 
variable only via two physically relevant functions: the strain of the solid and 
a properly normalized fluid mass density CiriUp et al. 20111 ICirillo et al. 20091 
ICirillo et al."20T0l [SciaFra et al. 2008] . "~ 

More precisely, consider the Jacobian x'si^Sji) of the solid placement map, 
which measures the ratio between current and reference volumes of the solid 
component, and let 

£(Xs,0:=[(Xs(^s,t))'-l]/2 (20) 

be the strain field. Let ^ : ^ R he the fluid reference density] we define 
the fluid mass density field 

mf(Xs, t) := g„j{^{Xs, t))^'{Xs, t) (21) 

Assuming that the mass is conserved, it is proven [Sciarra et al. 2008] that the 
field m£ can be interpreted as the fluid mass density measured with respect to 
the solid reference volume. 

We assume the total potential energy density <I> to depend on the fields 
and e and on their space derivative mj- and e' . Since = gQ^{4>)4>', m'^ = 



9 



+ Qoli^)^"^ £ = ((Xs)' - l)/2, and e' = XsXl we have that the 
Lagrangian density ^ = — $ depends on the space and time variables through 
the fields (/>, 4>' , (j)" , x's: and Xs- We then have 



SA 



dt 



( 



+ 



By using the above expression for the Lagrangian density we have that 

r-t; 



SA 



ti 



dt I dXs 



xdm^ d(j) dm'^ dcf) 

\dm^ d(j)' dm'^ d(j)' 

de 9$ de' 

V de dx's de' dx's 



dm'r 

! + — f 

dm'^ d(j)" 

, , a$ de' , „ 



(22) 



To get rid of the terms depending on the variations of the derivatives of the 
basic fields xs and (f) we integrate by parts. Via a standard computation we find 



5A= - 



dt 



dXs 



a$ dm^ d^ dm'^ 
dm^ d(t) dm'^ d(j) 

9$ dm^ a$ dm'f 



^drrif dcf. 
/9$ de 



(9$ dm'^ 



t2 



dt 



{( 



V de dx's de' dx's 
a$ dm'r Y2 /a$ de' ,N^2 
d^W ) ^d^dZf^^) 



dm'f d(f>' V dm'^ dcf)" 
de' 

Vfh' 



'f 

9$ de' 



de' dx's 



f^Xs} 



5$ 9m 



f 



dm 



5<j) 9toJ. / dm'f 



f 



(?£ 

V de dxs 



dm'^ dcj)' V dm'^ d<p' 
: a$ 9£' /d<P de' \'\ 



(23) 



^.5. Equations of motion 

By using the variational principle (|7| , the equation (16) for the virtual work of 
the drag forces, and the equation ( 23 1 for the variation of the action, we get the 
equations of motion 



Xs( 



de \de'J 



and £(0) 



r d<i> 


/ (9$ y 


-dm^ 


\dm'^/ . 



= R~Q' (24) 



and the boundary conditions 



I dm'r 



9£ 



(25) 
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We remark that, though they are partially written in terms of the fields 
and e, the equations (24) are evolution equations for the two fields xs and (j). 



Only under a brute approximation, see the geometrical linearization discussed 
in Section [Sj they will reduce to a set of evolutionary equations for the fields 
and £. 

Notice that those equations are not simply the one dimensional particulariza- 
tion of the equations (30)-(32) in the paper [Sciarra et al. 2008] . since there the 
reference mass density ^ was assumed to be constant. When infinite systems 
are considered not constant reference densities can be useful in order to have 
finite total fluid mass integral. In the sequel we shall always consider finite sys- 
tems so that it will be natural to assume the fluid reference density ^ : i?£ — > M 
to be constant. Under this hypothesis we get the following simplified expression 



Ik 



ide') 



and 



^?o.f 



r a$ 


(f.)'l 







R-Q' (26) 



for the equations of motion and 



9$ , . , /6 



S(t> 



a,7x^'5x^ + (|^~(|5)')xs^xs}^^ 



(27) 



= 



for the boundary conditions. We stress that from now on the symbol ^ denotes 
a real constant. 



3. Allen Cahn and Cahn Hilliard— like equations 



In this section we rewrite the equations of motion of the poromechanical system 
under the so called geometrical linearization assumption, namely, in the case 
when only small deformations are present in the system. We first introduce the 
displacement fields u{Xs,t) and 'w{Xs,t) by setting 



Xs{Xs,t)=Xs + uiXs,t) and (j){Xs,t) = Xs + w{Xs,t) 



(28) 



for any Xg e and i e M. We then assume that u and w are small, together 
with their space and time derivatives, and write the equations of motion up to 



the first order in u, w, and derivatives. By using (201, (21|, (17), and (18) we 
get 



^f = ^0 f(l + ^'); m := m,^- g^^= g^ fw' 



R sa Dw, and Q 



Sw' 
(29) 

where ~ means that all the terms of order larger than one have been neglected. 

We have introduced above the field m. In the following we shall imagine 
<i> as a function of m and m' and the equations of motion and the boundary 
conditions will be written in terms of this field. This can be done easily; indeed, 
since £ is constant, we have that d/dm^ = d/dm and = m! . From (26 1, 



(27), and (29), we get the equations of motion 



(1 



/9$ /a$ 



V 



and 



^^o,f 



9$ 
dm 



/ 9$ y 



Dw-Sw" (30) 
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and the associated boundary conditions 



f 9$ , d<^ ^ 

< — :oe+— — :dm- 



ide' 



dm' 



V dm V dm' 



S 

£»QfH m 



(31) 



A set of boundary conditions implying that (31 ) are satisfied is 



\di' 



Se 



d^ 
dm' 



9$ 

dm 



/ 9$ 
\dm' 



S 



d^ 

Ik 



= 



(32) 

where the notation above means that the functions in brackets are evaluated 
both in £i and £2- With this choice it is possible to fix the boundary conditions 
directly on fields m and e (and derivatives) . The equations ( 30 1 are evolution 



equations for the original kinematic fields, but, by treating separately the Darcy 
and the Stokes effects, we can deduce a system of equations for the fields m and 
e. 



The equations (31 1 are an extension to the case of second gradient elastic- 



ity of classical natural and essential boundary conditions [dell'Isola et al. 2009| . 
The first equation ( 32 1 is the additional boundary condition due to the presence 
of the gradient terms in the potential energy density This equation specifies 
essential boundary conditions on the derivatives of the displacement fields or 
natural boundary conditions on the so called double forces, see [Germain 1973| . 
The generalized essential boundary conditions can be read as a prescription on 
the derivative of the independent fields xs and 0, see equations (20) and (21 1; 



whilst the extended natural boundary conditions prescribe, on one hand, the 
additional forces which the solid continuum is able to balance at the bound- 
ary and, on the other, the wetting properties of the fluid which fills the pores 
[Seppecher 1989] , 



The second and the third equations ( 32 1 provide natural boundary condi 



tions prescribing the chemical potential of the fluid and the traction exerted 
on the overall porous solid, respectively. These conditions extend the classical 
ones (which can be easily found in the literature) for a porous solid suffering 
internal stresses due to applied tractions and for a saturating fluid with fixed 
chemical potential; see [Baek et al. 2004) for more details. In this case we use 
the words generalized traction and chemical potential because of the additional 
contribution to stress (d^/de) and chemical potential {d^/dm) provided by 
the spatial derivatives of the corresponding hyperstress fields, say {d^/de')' 
and {d^/dm'y. 

Pure Darcy dissipation. We consider the system of equations (301 for S — 0. 
Recalling that m = ^w' , by deriving the second between the equations of 
motion with respect to Xg, we get 



(1 



a$ /d<^\' 
^ \d^) 



= and 



^o,f 



\dm 



\dm' 



= Dm (33) 



Those equations, together with the boundary conditions (32), are a partial dif 



ferential equation problem for the two fields m and e. Note that, by exploiting 
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the third boundary condition ( 32 ) in £i, we get the PDE problem 

de 
/5$ 



\de' 



6e 



9$ 



<9$ 


/ 9$ y 




-dm 







9m' 



■Sm 



9m 



/ 9$ 
\dm' 



(34) 







Pure Stokes dissipation. We consider the system of equations (30) for D = Q. 



RecaUing that m = ^w' , the equations of motion (30) become 

-9$ 



/9$ /9$\' 
V9e ^ V9e'/ 



= and 



dm 



9$ 
9m' 



o,f 



= (35) 



Those equations, together with the boundary conditions (32), arc a partial dif- 
ferential equation problem for the two fields m and s. Note that, by exploiting 

H in. 

9$ 



the second and the third boundary condition ( 32 ) in , we get the PDE problem 

9$ /9$\' „ , 9$ 

dm 

9$ 



de 
/9$ 



V9e' 



6e 



dm 



= and 
6m 



/ 9$ Y 



5 
o.i 



m = 



(36) 



= 



It is very interesting to note that, provided the potential energy density $ 
depends on the space derivatives m' and e' of the fields m and e as a quadratic 
form with constant coefficients, it follows that the second between the equa- 
tions of motion (34) becomes a Cahn-Hilliard-like equation with driving field 
still depending parametrically on e, while the second between the equations of 
motion ( 36 ) becomes an Allen-Cahn-like equation with driving field depending 
parametrically on e. More precisely we specialize the model we are studying by 
choosing the second gradient part of the dimcnsionless potential energy, that is 
we assume 

$(m',£',m,e) ^[fci(e')^ + 2fc2£'m' + fc3(m')^] + *(m,e) (37) 

with fci, fcs > 0, fc2 G Hi such that kik^ — A:| > 0. These parameters provide 
energy penalties for the formation of interfaces; they have the physical dimen- 
sions of squared lengths and, according with the above mentioned conditions, 
provide a well-grounded identification of the intrinsic characteristic lengths of 
the one-dimensional solid and the fluid. 

With this choice of the second gradient potential energy, in the pure Darcy 
case the PDE problem ( 34 ) becomes 

(fcie" -I- k2m") = and Dm = -g^ H 



de 



(j(kie' + k2m')Se + {k2e' + kj,m')5r 



k2£" + fcsm" 
{k2e" + k^m" - 



9*\" 
dm J 



9* 
dm J 



while in the pure Stokes case the PDE problem ( 36 ) reads 
d-^ .5 



de 



[kie" + k2m") = and 



k2£" + fcsm," 



o,f 



9* 

dm 



Ukie' + k2m')de + {k2e' + kzm')5m\ = 

^ ' tl,t2 



■■ 
(38) 



(39) 
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We discuss, now, some general properties of the PDE problems (38) and 
(39). First of all we recall that an harmonic function of a single real variable 
equal to zero at the extremes of an interval is necessarily equal to zero in the 
whole interval. We then have that for both the two problems ( 38 ) and ( 39 ) the 
stationary solutions are the solutions of the stationary problem 

— - (/tie" + fcsm") = and fcse" + hm" ~^^0 
oe dm 



((kie' + k2m')Se + (k2e' + k-im')5m\ = 



(40) 



widely studied in [Cirillo et al. 2011[|Cmllo et al. 20091 ICirillo et al. 2010] . 

The dissipative character of the physical problem we are studying reflects in 
the existence of the not increasing energy functional 



<l)(TO',e',TO,£)dXs 



-[ki{e'Y + 2fc2£'m' + ks.{m'Y] + *(to, e) 



(41) 



dXs 



for both the systems ( 38 ) and ( 39 ) . To prove this we compute the time derivative 
of the functional F evaluated on the fields xsi^Sit) and (/)(Xs,t). We first get 



dF 
d^ 



{kie' + k2m')e' + {k2e' + fc3m')m' + - — m + — -e 

am oe . 



and, by integrating by parts, we obtain 

\ikie' + k2m')e + {k2£' + fc3rn')m]^^ 

- {kie" + k2m")e - (fcae" + k3m")m + ^—m 

dm 



dt 



(42) 

Assume, first, that the fields m and e are solution of the Allen-Cahn-like 
evolution equations (39), i.e., the PDE problems describing the evolution of the 
system when Stokes dissipation is considered. Moreover, assume that the bound- 
ary conditions have been chosen on the fields or on their first space derivatives 
so that those in ( 39 ) are fulfilled, we have that 

TfYdXs 



dF 
~dt 



(43) 



which proves that the functional F in not increasing along the motions of the 
Allen-Cahn-like system ( 39 ) . 

Assume, now, that the fields m and e are solution of the Cahn-Hilliard- 
like evolution equations (38), i.e., the PDE problems describing the evolution 
of the system when Darcy dissipation is considered. Moreover, assume that 
the boundary conditions have been chosen on the fields or on their first space 
derivatives so that those in (38) are fulfilled, from (42) we have that 



dF 
~dt 



D 



(fcae" + ksm") + 



dm 



{k2e" + Aigm") + 



dm 



dX. 
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By integrating by parts and exploiting the second boundary condition in ( 38 ) 
we get 



dF 

dT 



D 



{[-(fc2e" + A;3m") + |^]'} 



dXc 



(44) 



which proves that the functional F in not increasing along the motions of the 
Cahn-Hilliard-like system ( 38 ) . 



In conclusion, in this section, we have written the equations of motions ( 26 ) 
and the associated boundary conditions ( 27 ) under the geometrical linearization 
assumption. We have proven that, provided the Darcy and the Stokes effects 
are treated separately, those equations reduce respectively to the Cahn-Hilliard- 
like (38 1 and to the Allen-Cahn-like (39) PDE problems. Moreover, we have 
remarked that the stationary solutions of those problems are given by the same 
system of equations (40). Finally, we have proven that the energy functional 
( 41 ) does not increase along the solutions of both the two systems ( 38 ) and ( 39 1 , 
which suggests that in both cases the motions will tend asymptotically to the 
stationary profiles. 



4. Pressure driven phase transition 

We apply, now, the theory developed above to the special model introduced 
in ICirillo et al. 2011) ICirillo et al. 2009| ICirillo et al. 2010] and whose station- 
ary behavior has been widely discussed in those papers. We consider the fol- 
lowing expression for the total potential energy density in the perspective of 
describing the transition between a fluid-poor and a fluid-rich phase 

*(m,£) := ^m^(3m^- 86£m -f 66^e^) -I- *3(m,£) (45) 

where 

*g(m,e) := pe -I- -e^ -I- -a(m - fe)^ (46) 

is the Biot potential energy density |Biot 194T] . o > is the ratio between the 
fluid and the solid rigidity, 6 > is a coupling between the fluid and the solid 
component, p > is the external pressure, and a > is a material parameter 
responsible for the showing up of an additional equilibrium. 

In the papers jCirillo et al. 20111 ICirillo et al. 2009| ICirillo et al. 2010j we 



have studied the stationary solutions of the equations ( 38 ) and ( 39 ) , that is to 



say we have studied the problem ( 40 1 . In this section we give a brief account 
of those results and then wc will study numerically the not stationary solutions 
and, in particular, discuss how the stationary ones are approached. The two 
cases, pure Darcy and pure Stokes dissipation, will be discussed separately. 
First of all we write explicitly the stationary problem corresponding to the 



potential energy (45). By (40) we get 



-{2/3)abm^ + ab'^m'^e + p + e - ab{m~ be) - kis" - k2m" = 
am? — 2abm?e -I- ab'^rae^ + aim — be) — k2e" — k^m" = ^-^^-^ 

({kie' + k2m')Se + {k2e' + k^m')5ni\ = 
where the last line is the boundary condition. 
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4-1- Phases: constant stationary solutions 

In [Cirillo et al. 20lTllCirillo et al. 2009| we have studied the constant solutions 
of ( 47 ) which are called phases of the system. We have proven that there exists 
a pressure pc, called critical pressure, such that for any p G [0,pc) there exists 
a single phase (ms(p), es{p)), called the standard phase, which is very similar to 
the usual solution of the Biot model. For p > pc a second phase {m^{p),e^{p)), 
richer in fluid with respect to the standard phase and hence called fluid-rich 
phase, appears. 

We have shown that the standard phase (ms(p), es(p)) is the solution of the 
two equations m — be andp = fi{e), for any p > 0, where /1(e) — e— ab^e^/S. 
On the other hand the fluid-rich phase {■m^{p),e^{p)) is the solution, with the 
smallest value of e (recall e G (—1/2,0), so that the smallest value has indeed 
largest modulus) , of the two equations m = m+ (e) and p = f+ (e) , where 



4a 



and 



f+{e) :=— £ + a6[m+(£) — he] — a}p'em\[e) + -abm^j^{e) 



For e < -~2/{by^ a/a) the function /+(e) is positive, diverging to +00 for e 
—00, and has a minimum at £c such that /+(ec) = Pc] this explains why the 
fluid-rich phase is seen only for p > pc- 

Moreover it has been shown that for any p > the point (ms(p), es(p)) is 
a minimum of the two variable potential energy ^'(m,e) with p fixed, while 
{m^{p), £f(p)) is a minimum for p > pc and a saddle point for p = pc- For more 
details we refer to [Cirillo et al. 20lT] . 

4-2. Profiles: not constant stationary solutions 

In [Cirillo et al. 2010] it has been proven that there exists a unique value pco 
of the pressure, called coexistence pressure, such that the potential energy of 
the two phases is equal. More precisely, it has been proven that the equation 
'^{ms{p),es{p)) = \E'(mf(p), ef(p)) has the single solution pco- 

The behavior of the system at the coexistence pressure is particularly in- 
teresting; from now on we shall always consider p — pco and, for this reason, 
we shall drop p from the notation. When the external pressure is equal to pcoj 
none of the two above phases is favored and we ask if profiles connecting one 
phase to the other exist. More precisely, in [Cirillo et al. 2010] we have ad- 
dressed the problem of the existence of a connection between the two phases. 



that is, a solution of the stationary problem (471 on M tending to the standard 
phase for Xg — )■ —00 and to the fluid-rich one for Xg — )■ -l-oo. Using results in 
[Alikakos et al. 2008] we have proven that such a connection does exist provided 
kiks - /c| > 0. 

We have also shown that, for fcifca — fc^ — (degenerate case), the problem of 
finding a solution of the stationary problem can be reduced to the computation 
of a definite integral. Indeed, in such a case one performs the rotation of the 
Cartesian reference system 

m + ke , ~km + e 

X := , and y := — , (48) 
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in the plane m-e, where k :~ ^2/^3 = it-^/fci/fca, and defines 

V{x,y) = -^{m{x,y),e{x,y)) (49) 
Then one shows that the two fields m{Xs) and £{Xs) are solutions of the two 



equations (471 if and only if the corresponding fields x{Xs) and y{Xs) satisfy 

dV dV 
k3il + e)x" = - — ix,y) and — (x,2/) = (50) 

The root locus of the constraint curve dV{x,y)/dy = is made of a certain 
number of maximal components such that each of them is the graph of a function 
X e M — > y{x) € M; for each of them the first between the two equations 



( 50 1 becomes a standard one dimensional conservative mechanical system with 
Xg, interpreted as time, kinetic energy (1 + A;^)fc3(x')^/2, and potential energy 
V{x,y{x)). 

Since the function V has been obtained by flipping the sign of the function 
^' and rotating the coordinate axes, then at the coexistence pressure it has the 
two absolute maximum points {xs,ys) and (xj, yj) corresponding, respectively, 
to the standard and to the fluid-rich phases. Since (tos,£s) and (m^, e^) satisfy 
the equations 5'm(m,e) = and ^'^(to,^) = 0, we have that the two points 
(2^Sj2/s) and (x^,y^) are solutions of the constraint equation dV{x,y)/dy = 
and hence they belong to the constraint curve. 

In [Cirillo et al. 2010] we have seen that there exist values of the second 
gradient parameters fci, fc2j and ^3 such that the two points above fall on the 
same maximal component of the constraint equation. Since, in this case, the 
function V has two isolated absolute maximum points which, by hypothesis, 
belong to the same maximal component of the constraint curve, we have that 
the function V{x,y{x)) of x has two absolute isolated maxima in and x^. 
The motion of the equivalent one dimensional conservative mechanical system 
corresponding to the energy level l^nax '■= l^(a^Sj2/s) yields the heteroclinic 
orbit, that is to say it yields the above mentioned connection on M between the 
standard and the fluid-rich phases. 



With a similar approach one can find solutions of the system (47) with 
Dirichlet boundary conditions on a finite interval, say — Q and = 1- For 
instance, assume xs < x^ and consider the boundary condition a;(0) — xs and 
a::(l) — x^. Motions of the equivalent mechanical system started in xs at time 
with positive velocity a;'(0) = u > (note that if it were > x^ one should 
consider a negative velocity at time 0) will reach Xj? in a finite time, which is a 
decreasing function of the initial velocity. It is then possible to choose properly 
the initial velocity so that 




. , fc3(l + fc2) ^ 

where E — (l + A:^)fc3f;^/2 + Vrnax is the total energy of the motion of the equiva- 
lent one-dimensional conservative system. Once E has been found, the solution 



of the stationary problem (47) with Dirichlet boundary conditions m{Q) = mg, 
e(0) = £sj nn{l) = toj, and e(l) = is implicitly given, in terms of the variable 
X, by the definite integral 

_ ' fc3(l + fc^) 

' '^2[E-V{x,y{x))] 
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The solutions found in such a way for ki = 10~^, = 10~^, k-z = ±10"'^, 
a = 0.5, b ~ 1, a ^ 100, and at coexistence pressure, have been depicted in 
figures [T] and [2j 

As we have noted above, in the not degenerate case, i.e., kik^ — fc| > 0, in 
[Cirillo et al. 2010] we could prove the existence of a connection, but we were not 
able to find it explicitly. Studying in more detail the behavior of the solutions 



of the system (47 1 is a difficult task since, in this case, the problem cannot be 
rewritten as a one dimensional conservative mechanical system. We shall then 
study the behavior of the solutions close to the two phases by linearizing the 
stationary equations, while far from the points representing the phases we shall 
use a finite difference numerical approach. 

The first question that naturally arises looking at the solutions obtained in 
the degenerate case, see figure [l| is the nature of the bump. In the degenerate 
case the bump is due to the fact that the solution of the stationary problem, 
once projected onto the m-e plane, has to lie on the constraint curve, so that 
its shape depends on that of the constraint itself. In the not degenerate case no 
constraint curve exists, in other words it does not exist any a priori privileged 
manifold on the plane m-e. It is then possible (in principle) that the bump is 
associated to solutions oscillating close to the phases. 

In order to discuss the existence of such oscillatory solutions we linearize the 
problem close to the phases. Let {fh,e) be one of the two points (TOs,es) and 
(m£, ej) representing, respectively, the standard and the fluid-rich phase. By 



expanding the first gradient part of the stationary equations (401 around (m, e) 
to the first order, we get the linearized equations 

kie" + k2m" — ^'^^(m, e){m — rfi) — ^^^(m, e){e — e) — Q 
k2e" + k^m" - ^'™„(m,e)(TO - fh) - ^^^{fh,£){£ - e) = 

If we let qi = e — e and q2 — m — fh, the above equations become 

Tq" - *g = (51) 
where we consider the column vector g e and set 

^ and *='^ *e£(™,£) *£m(m,£) 



h ks J \ *TO£(m, e) '^rnm{'m,s) 

Note that the matrix T is trivially real and symmetric; moreover, in the not 
degenerate case fcifca — > it is also positive definite. The matrix >&, on 
the other hand, is real and, by Schwartz theorem, symmetric. Moreover, as 
we proved in | Cirillo et al. 2011] . it is positive definite at any pressure greater 
than the critical one pc] in particular it is positive definite at the coexistence 
pressure. 



Looking for solutions of (51 1 in the form q = cexp{AXs}, with c e 
and A G C, we get that the constant vector c has to be a solution of (A'^T — 
'4')c = 0. Set /i = A'^; in order to have not vanishing c solving the above 
equation, fi has to be a solution of the secular equation det(/xT — "if) — 0. 
Since T is real, symmetric, and positive definite and '4' is real and symmetric, it 
follows that the solutions fii and fi2 of the secular equations are real. Moreover, 
since ^ is positive definite we have that those solutions are positive themselves 
(see, for instance, [Gantmacher 19751 Section 40, Chapter 6]). Thus the secular 
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equations have four real pairwise opposite eigenvalues. This implies that the two 
fixed points are unstable, as it also follows from the existence of the connection, 
and that no oscillation around the fixed points is expected. 



We come, finally, to the numerical study of the stationary problem (47) in 
the not degenerate case on a finite interval (we shall always use the interval 
[0, 1]). We use the finite difference method powered with the Newton-Raphson 
algorithm. The substitution rules we adopted are standard and have been re- 



ported in appendix Appendix B We solved the stationary problem (47) with 
the Dirichlet boundary conditions rn(0) = mg, e(0) = £s, = Wf, e(l) = Ef, 
at the coexistence pressure for a — 0.5, b — 1, and a = 100. For the second 
gradient coefficient we have considered two sets; the first set is fci — 10~^, 
ks = 10^3, and fc2 = -0.4 x lO'^, 0.2 x IQ-^, 0.8 x lO'^ and the related numeri- 
cal results have been depicted in the bottom row of figure [T] and in figure [2] The 
second set is ki = 10"^, fcg = 10-^ and fca = -0.4 x lO^^, 0.2 x lO^^, 0.8 x 10"^ 
and the related numerical results have been reported in the top row of figure [T] 

On the bottom (resp. top) row of figure [l] we have depicted the m and 
e-profiles as functions of the space variable Xs; solid lines refer to the degen- 
erate case k2 = ±10~^ (resp. k2 = ±10~^), while dotted lines refer to the 
not degenerate cases ^2 = -0.4 x 10~^,0.2 x 10~^,0.8 x 10~^ (resp. fc2 = 
—0.4 X 10~^,0.2 X 10~^,0.8 X 10^^). In figure [2] the same solutions have been 
depicted on the plane m~e. The behavior of the solutions of the stationary 
problem in the degenerate and in the not degenerate case is similar. The exis- 
tence of the bump in the m and e-profiles is essentially due to the shape of the 
manifold of the m-e plane on which the solutions lie. 

It is worth noting that, even in the degenerate case, due to the asymmetry 
of the potential energy V, the position of the interface in the stationary profile 
depends on the value of the second gradient constants fci, k2, and k^. However 
in the limit for small values of ki, k2, and fcs it is observed that the interface 
tends to a definite position [Cirillo et al. 2012] . 

4-.3. Evolution under pure Stokes dissipation: Allen- Cahn-like equations 

The equations of motion describing the evolution of the system subjected to a 



pure Darcy drag force are the equations ( 39 ) which, for the special model ( 45 1 , 
read 

J ~{2/3)abm^ + ab'^m'^e + p + e - ab{m - be) - kie" - k2m" =Q 

I f I*^"^^ ~ 2a6m^£ + ab'^me^ + a{m — be) — k2e" — k^m"] = Sra ^ ' 

with the boundary conditions 

{{k2m! + kie')5e + {ksm' + k2e')Sm)i,^^i^ = (53) 

We have studied the above PDE problem with the same numerical approach 
used for the stationary problem and discussed in Section |4.2| That is we have 
used a finite difference method powered with the Newton-Raphson algorithm 
and with the standard substitution rules reported in [Appendix B[ 



In the figures [3||5| we have depicted the solutions of the problem ( 52 ) with 
Dirichlet boundary conditions m(0) = ms, e{0) — es, m(l) — m^, and e(l) = e^ 
on the finite interval [0,1], at the coexistence pressure for a = 0.5, 6=1, 
a = 100, fci = ^2 = /ca = 10~^. All the graphs refer to the same spatially 
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random initial condition on the interval [0, 1]. Details on the depicted times can 
be found in the caption of the related figures. 

Figures [3] and |4] show how the limiting stationary profile is formed. At 
times of order one a deviation from the initial random profile is noticed and 
a bumped profile with roughness due to the randomness of the initial con- 
dition is formed. The profile is smeared out in a time of order 10^ and an 
interface, far from the position of the stationary profile, is seen. The inter- 
face then undergoes a metastable, purely dissipative, evolution approaching 
the stationary state in times greater than 10^. This coarsening process has 
been widely studied for the standard Allen-Cahn equation, see for instance 
[Carr et al. 19891 IFusco et al. 19891 IWard 1994) . 

Figure [5] shows that the solution is soon attracted by the manifold in the 
plane m-e defined by the stationary solution. Afterwards the projection of the 
profiles onto the m-e plane keeps close to such a curve. We recall that, when 
the parameters ki, ^2, and k^, as in the cases illustrated in the figures [3]^ are 
chosen so that kik^ — A;| = 0, the algebraic equation of the limiting manifold is 



the second above the two equations ( 50 1 



We have measured the way in which the stationary solution is approached 



considering the difference between the value of the energy functional (41) at 
time t and that of the corresponding stationary solution. Results are reported 
in figure [6j We note that the energy functional is a decreasing function of time 
as proved in Section |3] (see equation (43)). 



Some of the features described above are peculiar of the initial random con- 
dition. We have also studied the same Allen-Cahn problem, with the same 
parameters and the same boundary conditions, with a deterministic stratified 
initial state. Results for the e-profile are depicted at times in the figure [7) we 
have not reported the m-profiles and the m-e-plane view since they do not add 
any more information. 

At times of order 10~^ a deviation from the initial stratified state is detected. 
The profile deforms into a new one with a droplet-like shape, at times of order 
10^, finally an interface profile nucleates at times of order 10^, as in the previous 
case corresponding to spatially random initial data. Similarly to this case the 
stationary interface is approached with the same velocity. Even in this case we 
have checked the way in which the solution tends to the stationary state via the 
energy functional (41). Results are reported in figure [sj 



4-.4- Evolution under pure Darcy dissipation: Cahn-Hilliard-like equations 
The equations of motion describing the evolution of the system subjected to a 



pure Darcy drag force are the equations ( 38 ) which, for the special model ( 45 ) , 
read 

{-~{2/3)abm^ + ab'^m'^e + p + e ~ ab{Tn — be) — kie" — k2m" = 
g'^ ^{am^ — 2abm'^e + ab'^me'^ + a{iTi — be) — k2e" — k^m")" — Dm ^ ' 

with the boundary conditions 



{{k2m' + kie')6e + {k^m' + k2e')Sm)e^j.^ — 

{am^ — 2abm^e + ab^me^ + a{m — be) — k2e" — k'im")t^^i.^ = 



(55) 



We have studied the above PDE problem with the same numerical approach 
used for the stationary problem and for the AUen-Cahn-like system. We have 
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repeated the same sequence of numerical computations, that is we have used 
the same parameters and the same initial and boundary conditions as for the 
above Allen-Cahn-like case. 

In the figures |9 11 we have depicted the solutions of the problem ( [54| with 



Dirichlet boundary conditions m(0) = mg, e(0) = Eq, m(l) = m^, and e(l) = 
on the finite interval [0,1], at the coexistence pressure for a = 0.5, 6=1, 
a = 100, ki ^ k2 — = 10~^. All the graphs refer to the same initial condition 
chosen randomly on the interval [0, 1] as that considered for the solution of the 
Allen-Cahn-like equations. Details on the depicted times can be found in the 
caption of the related figures. 

Figures [9] and [TO] show that, starting from a random spatial distribution of 
strain and fluid mass, micro-structured paths appear at times of order 10~^. 
They progressively evolve towards a droplet-like profile which is formed at times 
of order 10~^. The droplet undergoes a metastable evolution at times of order 
10^ up to the formation of an interface. This coarsening process has been well 
studied for the standard Cahn-Hilliard equation, see for instance [Ward 19981 
lAlikakos et al. 19911 IBates et al.HMl IBates et al.T995] . The last part of the 
evolution is a the motion of the interface approaching the stationary state at 
times of order 10^. 

Figure [Tl] shows that even the solution of the Cahn-Hilliard-like system 



(54), as that of the Allen-Cahn-like system (52), is attracted immediately by 
the manifold in the plane m-e defined by the stationary solution and, afterwards, 
the projection of the profiles onto the m-e plane keeps close to such a manifold. 
Note that the profile temporarily leaves the stationary manifold at time 3.7 in 
correspondence of the shrinking of the droplet present in the e-profile. 

As for the Allen-Cahn-like case, also in the Cahn-Hilliard-like case we have 
studied the way in which the stationary state is approached by using the energy- 



like functional (41) (see figure 12). 

Finally we have also studied the evolution of the system starting from a 
deterministic stratified initial datum; also in this case we used the same initial 
condition as the one used for the ( 52 ) problem, see figure [t] Results have been 
plotted in figures [13] and [T4j 

4-.5. Discussion of the numerical results 

Comparing the profiles of the strain and fluid mass density due to a pure Stokes 
(Cahn-Hilliard evolution) and a pure Darcy (AUen-Cahn evolution) dissipation, 
for spatially random or stratified initial conditions, some remarks naturally arise 
concerning short-time and long-time behavior. 

Monitoring the evolution process described by the Cahn-Hilliard partial dif- 
ferential equations and considering spatially random initial data, short-time 
micro-scale oscillations can be detected at times of order 10~^; for these char- 
acteristic times the profile is smeared out and progressively deforms firstly into 
a metastable drop-like shape, which is approached at times of order 10~^, and 
lastly into an interface-like graph which is reached at times of order 10""^. Anal- 
ogously the short-time evolution process, described by the Allen-Cahn dif- 
ferential equations and stemming from the same randomly distributed initial 
condition, still exhibits a kind of intermediate feature and finally provides an 
interface-like graph which propagates towards the stationary state. However 
this last is approached at times several order of magnitude greater than those 
characteristic of the Cahn-Hilliard evolution process (in particular ^ 10"*). 
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The same difference is appreciated when considering stratified initial data; 
again the short-time behavior provides for both evolutions a deformation of the 
initial profile into a drop-like new one at different time scales. Looking at the 
profiles related to different characteristic times one can assess that the evolution 
process governed by the Allen-Cahn-like set of partial differential equations is 
affected by initial conditions during a characteristic time interval much larger 
than the corresponding one typical of the Cahn-Hilliard-like equations, see 
figures [3]: and [9]:, which correspond to the formation of the drop-like shape of 
the profile. 

It is worth to notice that both the Allen-Cahn and the Cahn-Hilliard evolu- 
tion processes separately exhibit a consistent behavior when the initial data is 
changed. In particular, starting from spatially random or stratified initial condi- 
tions, both equations provide, after a transitory regime, a quasi-static evolution 
of an interface-like profile which nucleates and propagates at the same time and 
with the same velocity. On the other hand the numerical simulations suggest 
that the quasi-static evolution process starts at different characteristic times 
and develops with different velocities when the Darcy and Stokes dissipation 
mechanisms are considered. 

As already mentioned these two mechanisms correspond respectively to the 
presence (Darcy & Cahn-Hilliard regime) or absence (Stokes & Allen-Cahn 
regime) of scale separation and, thus, to the cases in which the fiuid dissipative 
flow is parametrised either by the permeability of the solid or by the viscosity 
of the fluid. In these two situations the micro-structure of the corresponding 
porous media is quite different: materials with scale separation are constituted 
by solid grains which do not hinder the fluid flow strongly; the grains can be 
regarded as undeformable and impermeable, suffering an equivalent drag force 
due to the average fluid velocity. On the other hand materials with no scale 
separation are micro-porous and the fluid flow is much more hindered by the 
micro-structure of the skeleton. This could provide a possible explanation for 
the very high discrepancy between the characteristic times of evolution for the 
Cahn-Hilliard and the Allen-Cahn-like system of equation: when dealing with a 
Cahn-Hilliard dissipative process the time necessary to approach the stationary 
state is lower than that necessary to reach the same state through an Allen- 



Cahn evolution, see figures M M [12] and 14 



Comparing the distance between the current and the stationary strain energy 
for the Cahn-Hilliard and Allen-Cahn evolution process, three or two different 
regimes can be appreciated. These indeed correspond to the smearing out of 
the random oscillation (if any), to the formation of the intermediate drop-like 
solution, and tothe nucleation of the interface-like profile. 

Deeper investigations will be developed in the future in order to better ex- 
plain such phenomenon, also accounting for different boundary conditions. Par- 
ticular attention should be devoted to the study of initial boundary value prob- 
lems involving Neumann boundary conditions, which naturally allow to describe 
the behavior of granular media under consolidation. 
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Figure 1: Bottom: solutions {s{Xs) on the left and m{Xs) on the right) of the stationary 
problem ( |47[ l with the Dirichelet boundary conditions m(0) = mg, e(0) = Ss, m{l) = rrif, and 
e(l) = ££ on the finite interval [0, 1], at the coexistence pressure for a = 0.5, 6=1,0 = 100, 
fci = 10-3, ks = 10-3, k2 = ±10-3 (solid lines), and k2 = -0.4 X 10-3, 0.2 X IQ-^, 0.8 X IQ-^ 
(dotted lines). Top: the same for fci = 10— k^ = 10—^, ^2 = itlO-'^ (solid lines), and 
k2 = -0.4 X 10-2,0.2 X 10-2,0.8 X 10"^ (dotted lines). 



Appendix A. Stokes virtual dissipation in terms of the basic fields 



Derivation of the equation ( 15 1 for the virtual work of the dissipation due to 
the Stokes effect. We first note that from the equaUty (10) we have that 

(A.l) 

where in the first step we used the chain rule for the derivative of a composed 
function and in the second step the expression of the derivative of the inverse 
function. In the second term of the equation above the notation means that 
the derivative of (f>Xs/4'' is computed with respect to Xg, and the result is then 
evaluated in Xs = Xs^i^^^)- 

Now, we note that by evaluating the expression ^ for the variation of the 
field Xf in ^s = Xs^i^7 fo'" ^^^y x € Bt and i G M, we get 

SxiiXi\x,t),t) - 5xsiXs\^,t),t) = -^|^^|^^50(Xs (A.2) 
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Figure 2: Plot on the m-e plane of the solutions of the stationary problem j47| already 
depicted in the bottom row of the figure [l] 



As before, from (A.2) we get 



(A.3) 



variables x = xsiXs,t) we get (15) 



Finally, by inserting (A.l) and (A.3) in (14) and by performing the change of 



Appendix B. Finite difference approximations 

In this appendix we collect the finite difference substitution rules that we have 
adopted in our numerical computations. Let n a positive integer number. Let 
a — 1/n and t G M be respectively the space and the time increments. The 
space interval [0, 1] is subdivided into n small intervals of length a. Given a 
field h{Xs, t), for any i = 1, . . . , n — 1 and j G N, we set 



h'{ia,jT) « —[h{{i + 1)ct,jt) - h{{i - 1)(T, jr) 

ACT 



For the second space derivative we set 



1 



h"{ia,jT) « —[h{{i + l)a,jT) - 2h{iG,jT) + h{{i - l)a, jr)] 



for i = 1, . . . , 71 — 1 and j E N. For the third space derivative we set 

h"'iza,jT) « + 2)a,jT) - 2h{ii + l)a, jr) 

+ 2h{{i~l)<j,jT)-h{{i~2)a,jT)] 
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for i = 2, . . . , n — 2 and j e N. For the fourth space derivative we set 

h""{ia,jT)^^[h{{i + 2)a,jT) - Ah{{i + l)a,ir) + 6h{ia,jT) 

- 4h{{i - l)a, jr) + hiii - 2)a, jr)] 

for i = 2, . . . ,n — 2 and j G N. For the time derivative we finally set 
h{ia,jT) w - [h{ia,jT) - h{ia, {j - 1)t)] 

T 

for i = 0, . . . , n and j = 1, 
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Figure 3; Profiles (solid thin lines) e{Xs,t) obtained by solving the AUcn-Cahn-like system 
\52\ with a random initial condition and Dirichelet boundary conditions m(0) = mg, £{0) = esi 
m(l) = mj, and e(l) = ££ on the finite interval [0, 1], at the coexistence pressure for a = 0.5, 
b = I, a = 100, ki = k2 = = 10~^. The solid thick line is the corresponding stationary 
profile. Profiles at times t = 2, 5.4, 7, 15, 200, 50000 are depicted in lexicographic order. 
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Figure 6: For the same problem as in figure [s] the difference of the energy l |41| at time t 
and that corresponding to the stationary profile is reported as function of time. The vertical 
thin lines denote the times t = 2, 5.4, 7, 15, 200 at which the e and m-profiles are depicted in 
figures Islfs] 
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Figure 7: Profiles (solid thin lines) £(Xs,t) obtained by solving the Allen— Cahn— like system 
l |52| with a deterministic stratified initial state and Dirichelet boundary conditions m(0) = ms, 
e(0) = £s, m(l) = m£, and £(1) = ££ on the finite interval [0, 1], at the coexistence pressure 
for a = 0.5, b = 1, a = 100, fci = fc2 = fcs = 10~^. The solid thick line is the corresponding 
stationary profile. Profiles at times t = 0.07, 2, 10, 320, 340, 50000 are depicted in lexicographic 
order. 
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Figure 8: For the same problem as in figure [t] the difference of the energy l |41[ l at time t and 
that corresponding to the stationary profile is reported as function of time. The vertical thin 
lines denote the times t = 0.07, 2, 10, 320, 340 at which the e profiles are depicted in figures [t] 
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Figure 9: Profiles (solid thin lines) e{Xs,t) obtained by solving the Cahn-Hilliard-likc system 
l |54| with the same random initial state as the one used in figure [s] We used Dirichelet bound- 
ary conditions m(0) = ms, e(0) = es, m(l) = m£, and e(l) = ££ on the finite interval [0, 1], at 
the coexistence pressure for a = 0.5, b = 1, a = 100, fci = A:2 = fcs = 10~^. The solid thick 
line is the corresponding stationary profile. Profiles at times t = 0.009,0.04,0.2,3.7,3.9,89 
are depicted in lexicographic order. 
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Figure 12: For the same problem as in figure [9] the difference of the energy at time t 

and that corresponding to the stationary profile is reported as function of time. The vertical 
thin lines denote the times t = 0.009,0.04,0.2,3.7,3.9,89 at which the e and m-profilcs are 
depicted in figures [9j |ll| 
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Figure 13: Profiles (solid thin lines) e{Xs,t) obtained by solving the Cahn— Hilliard-like 
system ( |54| with a deterministic stratified initial state (the same used in figure [tJ and 
Dirichelet boundary conditions m(0) = ms, £(0) = es, m(l) = m£, and e(l) = ££ 
on the finite interval [0,1], at the coexistence pressure for a = 0.5, b = I, a = 100, 
ki = k2 = ks = 10~^. The solid thick line is the corresponding stationary profile. Pro- 
files at times t = 0.001, 0.01, 10.3, 11, 11.1, 89 are depicted in lexicographic order. 



38 



10- ■ 



Figure 14: For the same problem as in figure [T3] the difference of the energy (jJTj at time t and 
that corresponding to the stationary profile is reported as function of time. The vertical thin 
lines denote the times t = 0.001, 0.01, 10.3, 11, 11.1, 89 at which the e profiles are depicted in 
figures |13[ 
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